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Abstract. To solve numerically boundary value problems for parabolic 
equations with mixed derivatives, the construction of difference schemes 
with prescribed quality faces essential difficulties. In parabolic problems, 
some possibilities are associated with the transition to a new formulation 
of the problem, where the fluxes (derivatives with respect to a spatial 
direction) are treated as unknown quantities. In this case, the original 
problem is rewritten in the form of a boundary value problem for the sys- 
tem of equations in the fluxes. This work deals with studying schemes 
with weights for parabolic equations written in the flux coordinates. Un- 
conditionally stable flux locally one-dimensional schemes of the first and 
second order of approximation in time are constructed for parabolic equa- 
tions without mixed derivatives. A peculiarity of the system of equations 
written in flux variables for equations with mixed derivatives is that there 
do exist coupled terms with time derivatives. 



1 Introduction 

Investigating many applied problems, we can consider a second-order parabolic 
equation vifith mixed derivatives as the basic equation. An example is diffusion 
processes in anisotropic media. In desining various approximations for the corre- 
sponding boundary-value problems, we focus on the inheritance of the primary 
properties of the differential problem during the construction of the discrete 
problem. 

Locally one-dimensional difference schemes are obtained in a simple enough 
way for second-order parabolic equations without mixed derivatives [Tl[5] . Mixed 
derivatives complicate essentially the construction of unconditionally stable sche- 
mes of splitting with respect to the spatial variables for parabolic equations with 
variable derivatives, even for two-dimensional problems. 

In some problems, it is convenient to use the fluxes (derivatives with respect 
to a spatial direction) as unknow quantities. This idea may be implemented in the 
most simple manner for one-dimensional problems [3^. To introduce fluxes, mixed 
and hybrid finite elements are applied [US] . The original parabolic equation with 
mixed derivatives may be written as a system of equations for the fluxes. The 
basic peculiarity of this system is that the time derivatives for the fluxes in 
separate equations are interconnected to each other. For the problem in the flux 
variables, unconditionally stable schemes with weights are developed. Locally 
one-dimensional schemes are proposed for problems without mixed derivatives. 



2 Differential problem 

In a bounded domain J7, the unknown function u{x,t), x = {xi,X2, ■■■,Xm), 
satisfies the equation 



Q:,p — 1 ^ 

Assume that the coefficients k^js^ a, /3 = 1, 2, m satisfy the conditions 

mm m 
Q — 1 a,/3 — 1 CK— 1 

= ^/3q, a, /3 — 1, 2, m, X ^ f2 



(1) 



(2) 



for any ^0(0;), a — 1,2, m with constant k>0. Consider the boundary value 
problem for equation ([T]) with homogeneous Dirichlet boundary conditions 

u{x, t)=0, xedn, 0<t<T (3) 

and the initial conditions in the form 

u{x,0) =u°{x), xen. (4) 

We introduce a vector quantity q = (gi, 52, Qm)"^ (the index T denotes 
transposition) such that 

q — — /Cgradu, (5) 

where /C= (fca^) is a square matrix toxto (/C G K™™) with elements A:Q^(a;), a,/3 
1, 2, rn. Using this notation, equation ((T| may be written as 

Ou 

— +diYq = f,xef2,0<t<T. (6) 
ot 

We can write the above problem ©-(IS]) in the operator form. Scalar functions 
are considered in the Hilbert space H = L2{fi) with the scalar product and norm 
defined by the rules 



{u,v) — / u{x)v{x)dx, ||u|| = (u, u)"'"/^. 
For vector functions, we use the Hilbert space V = L2(J7), where 
(q,g) = y'/ 'laix)ga{x)dx, ||q|| (q,q)^^^ 

Taking into account we can treat the matrix /C as a linear, bounded, 
self-adjoint, and positive definite operator in V: 

/CiV^V, /C = /C*, k£<IC<k£, fc > 0, (7) 



where £ is the identity operator in V. Suppose Vu — — gradu,i.e. 

d d d ^ ^ 



V .v. ^V, V y , , . . . , j . (8) 

On the set of functions that satisfy the boundary conditions for the gradient 
and divergence operators, we have 

udrvqdx+ / qgradu da; = 0. 

It follows from this that 2?*q — divq, i.e., 

2?*:V^-H, P*^ fA,^,..., (9) 

\dxi 3X2 OXmJ 

In the above notation ([T])-®, from dS])-®, we obtain the Cauchy problem 
for the system of operator-differential equations 

^+V*q=fit), 0<t<T, (10) 

q = /CX>M, (11) 
u{0) = u°. (12) 
For the problem (IT])-(|4]), the following equation corresponds 

— + V*KVu = f{t), 0<t<T, (13) 

wich is supplemented by the initial condition (|12p . Taking into account that 

d ,^ ,^ d 
— /C-/C— =0, 
dt dt 

it is possible to eliminate u from the system of equations (fTO|) . (|lip that gives 
+VV*q = Vf, C = IC-\ 0<t<T. (14) 

In view of pTj) and (fT^ . we put 

q(0) = q° = ICVu°. (15) 

In constructing locally one-dimensional schemes (schemes based on splitting 
with respect to spatial directions), we focus on the coordinatewise formulation 
of equations HU]), (HI]), dH]) ) and ([H]). Let 

P = (2?1,I?2, • • ■ ,2?™)^, a; = (/Ca/3), C = (Ca/?), 



then the basic system of equations pH)) . ([TT]) takes the form 

^ + J2v:q^^f{t), 0<t<T, (16) 



Qa ^^ICa0T>i3U, a = l,2,...,m. (17) 

13=1 

The equation (fT3l) for u is reduced to 

^ + ^ i?;/c„/3p^w = /(t), < t < r. (18) 



For the flux components (see ([T4)) ). we obtain 

^Ca/3^ + ^I?„2?^<7^ 0<t<r. (19) 

/3=1 /3=1 

The equations of the system (|19p are connected with each other, and, more- 
over, the time derivatives are interconnected. The problem p2|) . (fTSl) seems to 
be much easier — we have a single equation instead of the system of m equa- 
tions. Nevertheless, some possibilities to design locally one-dimensional schemes 
for the system of equations are still there. 

Here we present elementary a priori estimates for the solution of the above 
Cauchy problems for operator-differential equations, which will serve us as a 
checkpoint in the study of discrete problems. Multiplying equation ([1^)) scalarly 
in 'H by u, we obtain 

\\u\\^^\\u\\ + ilCVu,Vu)^if,u). 
Taking into account (O and 

if,u)<\\f\\\\ul 

we arrive at 

|lHI<ll/l|. 

This inequality implies the estimate 

\\um<\\u°\\+ fwfmde (20) 

Jo 

for the solution of the problem [T2|) . (ITSt . 

Now we investigate the problem ([T4l) . ([TSl) . By the properties ([7]) of the 
operator /C, for C, we have 

C:V^V, C = C*, cS<C<c£, c = 'fc"^>0, c = fc"\ (21) 



In view of (ICT . we define the Hilbert space Vc, wfiere tlie scalar product and 
norm are 

(q,g)c = (Cq,g), ||q||c = (q,q)y^ 
Multiplying equation (ITSI) scalarly in V by q, we obtain 

l|q||c|||q||c + (2?*q,I?*q) = (2?/,q). 

In view of 

(P/,q) < ||2?/|k||q||c, 
we arrive at a priori estimate 

||q(t)||c<||q°||c+ f\\Vf{e)Ud0 (22) 
Jo 

for the solution of the problem (ITi)) , ([T5|) . 
3 Approximation in space 

We conduct a detailed analysis using a model two-dimensional parabolic problem 
in a rectangle 

n ~ {x \ X = {xi,X2), < Xa < la, « = 1, 2}. 

In ^2, we introduce a uniform rectangular grid 

Zj = {X \ X = {xi,X2), Xa=iaha, ia = 0, I, ... , Na, N^ha ^ la} 

and let uj be the set of interior nodes (ZU = cj U duj). On this grid, scalar grid 
functions are given. For grid functions y{x) = 0, a; G duj, we define the Hilbert 
space H — L2{uj) with the scalar product and norm 

(y, ^) = Y1 Vix)w{x)hih2, \\y\\ = {y, yfl"^. 

To determine vector grid functions, we have two main possibilities. The first 
approach deals with specifying vector functions on the same grid as it used for 
scalar functions. The second possibility, which is traditionally widely used, e.g., 
in computational fluid dynamics, is based on the grid arrangement, where each 
individual component of a vector quantity is referred to its own mesh. Here we 
restrict ourselves to the use of the same grid for all quantities, in particular, for 
setting the coefficients kapix), a,f3 = 1,2,..., m. 

Consider approximations for the differential operators 



We apply the standard index-free notation from the theory of difference schemes 
[B] for the difference operators: 

u{x + h) — u{x) u{x) — u{x — h) 

Ux — "J^ : U-x — • 

If we set the coefficients of the elHptic operator at the grid points, then 

LciaV ~ ~'^(kaaUxa)xa ~ '^i}^aaUxo,')xc,i OL — \,2. (23) 

More opportunities are available in approximation of operators with mixed 
derivatives. As the basic discretization [6], we emphasize 

Lalv = -\{kaf}UxJxff - ]^{kal3UxJxp, (24) 

L^lv = -^^{kapUxJxfi -]^[KfsUx^)xf,, a,/3 = l,2, a ^ /3. (25) 

Instead of L'^^p^-^'ap^ '"^'^ take their linear combination. In particular, it is 
possible [7] to put 

= + «,/3 = l,2, (26) 

In the general case, we set 

i";3 = x4'j + (1 - a,/3 = l,2, x = const . (27) 

The introduced discrete operators approximate the corresponding differential 
operators with the second order: 

L^o^u^ Uo.u + 0{hl), L^pu^ Ufi + Oih'), /3^a, a,/3 = l,2, (28) 

where — hf + h\ . 

We define a grid subset ZU, where the corresponding components of vector 
quantities are defined. Let 

bJ^ = {X \ Xi ^ il, ii = 0, 1, iVi - 1, X2= 12^12, «2 = 1, 2, ...,iV2 - 1}, 

= {a; I xi = «i, ii = 1, 2, A^i, X2 = 12/12, ^2 = 1, 2, N2 - 1}, 
io'^ = {x\ xi = iihi, ii = 1, 2, A^i - 1, X2 = 12/12, «2 = 0, 1, N2 - 1}, 
u}2 = {x \ xi ^ iihi, ii = 1, 2, A^i - 1, X2 = 12^12, ^2 1, 2, A^2}, 

and 

UJ — U!^ U UJi U UJ2 U UJ2 ■ 

For the grid vector variables, instead of two components, we will use four com- 
ponents, putting 

iqt,qi:Qt^Q2 f^ <lt^lt{x), a;ea;^, a = 1,2. 



For the grid functions defined on grids w^, a — 1,2, we define the Hilbert 
spaces H^, a — 1,2, where 

^ y(a;)^i;(a;)/^l/i2, \\y\\t ^ iiy,y)ty^\ a = 1,2. 

For the grid vector functions in F = © © H2 © H2 , we set 
2 

(q,g) = E((9"'ff^)^ + llq|l = (q,q)'/'- 

Now we construct the discrete analogs of differential operators Da, 2?^, a = 
1,2 introduced according to ([S]), (0). Using the above difference derivatives in 
space, we set 

Dtv = -V.^, x^ujI, a = 1,2, (29) 

so that D+ : H ^ H+, a = 1,2. Similarly, we define D' : H ^ H' , a = 1,2, 
where 

D^y^-ys^, xecu-, a = 1,2. (30) 

Thus 

D:H-^V, D = {D+,D^,D+,D2f. (31) 
For the adjoint operator, we have 

D*:V-^H, D* = {{D+r,{D^r,{D+r,{D2r), (32) 

and 

(D+r-.H+^H, {Dt)q = q,^, (33) 

[D-y-.H-^H, {D-)q^q^^, xetu, a = 1,2. (34) 

The above discrete operators approximate the corresponding differential opera- 
tors with the first order: 

Dtu:^V,,u + 0{h^), {DtTu = Vlu + 0{K), a = 1,2. (35) 

For the operator-differential equation (|13p . we put into the correspondence 
the equation 

^+D*KDy = ^it), 0<y<T, (36) 

where, e.g, (p{t) = f{x,t), x E uj. For equation (pG)) . we consider the Cauchy 
problem 

2/(0) = u". (37) 

The construction of the operator K is associated with the approximations 
(P5|) - (P7)) . The most important properties are self-adjointness and positive def- 
initeness of the operator K. The equation (|36p approximates the differential 
equation ([T3| with the second order. 



The system of equations ([TUl) . ([TT|) is attributed to the system 

^+i?*g = ^(i), 0<t<T, (38) 

g = KDy. (39) 
For the flux problem (1141) . (|15p. we put into the correspondence the problem 



C^+DD*g^ D(p{t), C = K-^, 0<t<T, (40) 

g(0) = KDu°. (41) 

Similarly to ([20)) . we prove the following estimate for the solution of the 
problem ([55 1) . ([57 ]) : 

||y(i)||<||^i"||+ fM0)\\d9. (42) 







For the estimate (|22]) , we put into the correspondence the estimate 



Wsmc < WDu'Wk + f \\D^{6)\\Kde (43) 
Jo 



for the solution of the problem pUl) , (|4T 



4 Operator-difference schemes 

We introduce a uniform grid in time with a step r and let = j/(i"'), = nr, 
n = 0, 1, -/V, A^T = T. For numerical solving the problem ([551). ([57|) . we apply 
the standard two-level scheme with weights, where equation is approximated 
by the scheme 

. n+l _ n 

^ ^ + A(ay"+i + (l-a)2/") = ^", n = 0, 1, ^ - 1, (44) 

r 

where 

A = D*KD, A A* > (45) 

and, e.g., ip^ — f{(Tt^^^ + (1 — cr)t"). Taking into account ([ST]), the operator- 
difference equation is supplemented with the initial condition 

2/° = (46) 

The truncation error of the difference scheme (|44)) - (l46l) is C'(|/i|^-|-r^-|-(CT— 0.5)t). 

The study of the difference scheme is conducted using the general theory of 
stability (well-posedness) for operator-difference schemes 0[8] . Let us formulate 
a typical result on stability of difference schemes with weights for an evolutionary 
equation of first order. 



Theorem 1. The scheme i44^ - (TS^ is unconditionally stable for a > 0.5, 
and the difference solution satisfies the levelwise estimate 

< ||y"||+r||<^"||, n = 0,l,...,^-l. (47) 
From ([Tr|) . in the standard way, we get the desired stabihty estimate 

\\y-+'\\<\\u^ + j2ry'i 

k=0 

which may be treated as a direct discrete analogue of the a priori estimate (|20p 
for the solution of the differential problem ([12]) , (|18l) . 

Schemes with weights for a system of semi-discrete equations ([551) . (15^ are 
constructed in a similar way. We put 

n+l _ n 

^ ^+i?*(ag"+i + (l-(7)g")-V'", n = 0,l,...,iV-l, (48) 

r 

g" = i^L»y", n = 0, l,...,iV. (49) 

The scheme (|48l) , (|48)) is equivalent to the scheme (|44)) . In view of Theorem 1 , it 
is stable under the restriction a > 0.5, and for the solution of difference problem 
(|i5|) . ([iSjl . the a priori estimate (gT]) holds. 

The special consideration should be given to the flux problem (PO]) . (PT|) . To 
solve it numerically, we apply the scheme 

n+l _ n 

^ + i?i5*((7g"+i + (l-a)g") = i?(^", n = 0,l,...,A^-l, (50) 

r 

g° = XDuO. (51) 

Theorem 2. T/ie difference scheme li5U\) . liSlfl is unconditionally stable for 
a > 0.5, and the difference solution satisfies the estimate 

||g"+ic<||g"||c + Tp(p"|k, n = 0,l,...,iV-l. (52) 
From (l52l) . it fohows the estimate 

n 

||g"+i||c < \\Du"\\k + J2 ^WDv'Wk, n = 0, 1, N-1, 

k=0 

which corresponds to the estimate for the solution of the problem (HUl) . (HT|) . 

The computational implementation of the unconditionally stable operator- 
difference schemes (P ^ " P5)) for the parabolic equation ([T|) with mixed deriva- 
tives is based on solving discrete elliptic problems at every time step. For the 
problem (j36p , (1371) , it seems more convenient to employ additive schemes (operator- 
splitting schemes) that provide the transition to a new time level using simpler 



problems associated with the inversion of the individual operators Z?*Z?q, a = 
1, 2 rather then their combinations. By the nature of the operators -D* , Da, a = 
1, 2, in this case, we speak of locally one-dimensional schemes. 

The issues of designing unconditionally stable locally one-dimensional schemes 
for a parabolic equation without mixed derivatives have been studied in detain. 
For parabolic equations with mixed derivatives, locally one-dimensional schemes 
were constructed in several papers (see, e.g., [9lfT0]V Strong results on uncondi- 
tional stability of operator-splitting schemes can be proved only in a uninterest- 
ing case with pairwise commutative operators (the equation with constant coef- 
ficients). For our problems the construction of locally one-dimensional 
schemes requires separate consideration. 

Let us investigate approaches to constructing locally one-dimensional schemes 
for the problem (HUl) . (HI]). The computational implementation of the scheme 
with weights ((5(I| . ([ST]) , which is unconditionally stable for a > 0.5, is associ- 
ated with solving the system of difference equations for four components of the 
vector g"+^. The equations of this system are strongly coupled to each other, 
and this interconnection does exist not only for the spatial derivatives (operators 
DfD2*, of of*), but also for the time derivatives (fci2 = fc2i 7^ 0). Thus, we 
need to resolve the problem of splitting for the operator at the time derivative, 
too. 

The simplest case is splitting of the spatial operator without coupling the 
time derivatives. Such a technique is directly applicable for the construction of 
locally one-dimensional schemes for parabolic equations without mixed deriva- 
tives, where 

kap{x) = kpa{x) = 0, a (3 — 1,2, ...,m, x E f2 (53) 

in equation ([T|). 

Assume that 

R = DD*, Q = dmg{D+{D+y,D^{D^y,D+{D+y,D^{D^Y), 

i.e., Q is the diagonal part of R. For numerical solving the problem PO)) . (HI]), 
we employ the difference scheme, where only the diagonal part of R is shifted to 
the upper time level. In our notation, we set 

n+l _ n 

— 5- + Q(ag"+i + (1 - a)g") + (i? - 0)g" = D^'\ n = 0, 1, iV - 1, 

(54) 

with the initial conditions according to (|5ip . 

Theorem 3. The difference scheme h51]) . is unconditionally stable for 
a > 2, and the difference solution satisfies the estimate 

l|g"+iB < ||g"||B+rp(p"||B-i, n = 0,l,...,iV-l, (55) 



where 



B^C + gtP - ^R. 



The scheme ([5T|) . ([51)1 has the first-order approximation in time. It seems 
more preferable, in terms of accuracy, to apply the scheme that is based on the 
triangular decomposition of the self-adjoint matrix operator R: 

R = Ri + R2, Rl = R2. (56) 

For the problem (|40p . (PT|). we construct the additive scheme with the splitting 
([5S)) . where 

71+1 _ n 

{C + <TTRi)C-^{C + crTR2)- — + Rs"" = Dip'' , n = Q,l,...,N-l. (57) 

T 

The main result is formulated in the following statement. 

Theorem 4. The difference scheme Ii51\] . 1561) - J??] ) is unconditionally stable 
for (7 > 0.5, and the difference solution satisfies the estimate h55]) with 

B ={C + aTRi)C^\C + (Tri?2) " ^R- 

The alternating triangle operator-difference scheme (|5T|) . (|56l) -(|57 |) belongs 
to the class of schemes that are based on a pseudo-time evolution process — the 
solution of the steady-state problem is obtained as a limit of this pseudo-time 
evolution. It has the second-order accuracy in time if cr = 0.5, and ony the first 
order for other values of a. 
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